Array measurements ROBAT 226.238:

• 5 mic array adafruit 5 i2s
• 48khz
• 1.5 meters distance (~ far field = 10𝛌; 𝛌max = fmin; 343/2 = 171 Hz --> 10𝛌 = 1715 Hz; at 1000Hz: 4.5𝛌 = 4.5*0.343 = 1.5435 m)
• measured clockwise (R: 0-180; L: 180-360)
• array rotated around the central mic kept a 1.5m from source

HW settings:

• Harman kandon AWR 445 vol = -25 db (front: 270 to 90)
• Harman kandon AWR 445 vol = -5 db (back: 100 to 260)
• fireface analog out 1/2 stereo vol = 0db
• tweeter #1
• Ref mic: gras +30 db fireface channel 9, +20db channel A power module

NOTE:

from IPython.display import Image
from IPython.display import display

display(
    Image(filename="./PXL_20250509_10061094.jpg", width=300),
    Image(filename="./PXL_20250509_10062195.jpg", width=300),
    Image(filename="./rme802_matrix.png", width=700),
    Image(filename="./rme802_mixer.png", width=700)
)

Output signal

Out signal used during the recordings: 1-40Khz 5ms sweep.

import scipy.signal as signal 
import numpy as np 
import matplotlib.pyplot as plt 
import sounddevice as sd
import soundfile as sf

# Set this variable to True to print, False to skip printing
do_print = False

#%%
# make a sweep
durns = np.array([3, 4, 5, 8, 10] )*1e-3
fs = 48000 # Hz

chirp = []
all_sweeps = []
for durn in durns:
    t = np.linspace(0, durn, int(fs*durn))
    start_f, end_f = 1e3, 20e3
    sweep = signal.chirp(t, start_f, t[-1], end_f)
    sweep *= signal.windows.tukey(sweep.size, 0.95)
    sweep *= 0.8
    sweep_padded = np.pad(sweep, pad_width=[int(fs*0.1)]*2, constant_values=[0,0])
    all_sweeps.append(sweep_padded)
    chirp.append(sweep)
    

sweeps_combined = np.concatenate(all_sweeps)

# Read the saved WAV file
out_signal, _ = sf.read('1_20k_5sweeps.wav')

# Plot the time-domain signal and spectrogram
plt.figure(figsize=(30, 8))

# Time-domain plot
plt.subplot(2, 1, 1)
plt.plot(np.linspace(0, len(out_signal) / fs, len(out_signal)), out_signal)
plt.title('Output Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')

# Spectrogram plot
plt.subplot(2, 1, 2, sharex=plt.gca())
plt.specgram(out_signal, Fs=fs, NFFT=512, noverlap=256, cmap='viridis')
plt.title('Spectrogram of the Output Signal')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.ylim(0, 25e3)
plt.tight_layout()

# plot the chirp
plt.figure(figsize=(30, 15))
for i, sweep in enumerate(chirp):
    plt.subplot(len(chirp), 2, 2 * i + 1)
    plt.plot(np.linspace(0, len(sweep) / fs, len(sweep)), sweep)
    plt.title(f'Sweep {i+1}')
    plt.xlabel('Time [s]')
    plt.ylabel('Amplitude')
    # Spectrogram plot

    plt.subplot(len(chirp), 2, 2 * i + 2)
    plt.specgram(sweep, Fs=fs, NFFT=32, noverlap=16, cmap='viridis')
    plt.title('Spectrogram of the Output Signal')
    plt.xlabel('Time [s]')
    plt.ylabel('Frequency [Hz]')
    plt.ylim(0, 25e3)

plt.tight_layout()
plt.show()

/home/alberto/miniconda3/envs/robat2/lib/python3.11/site-packages/matplotlib/axes/_axes.py:7947: RuntimeWarning: divide by zero encountered in log10
  Z = 10. * np.log10(spec)

Analysis

1. Angular audiofile import and extraction from multiwav to single channel recordings

# %% Libraries and files
import os
import soundfile as sf

# Load audio files, then plot a 6x6 grid
DIR = "./2025-05-09/original/"  # Directory containing the audio files
audio_files = os.listdir(DIR)  # List all files in the sweeps directory
audio_files.sort()  # Sort the files in ascending order

# Directory to save the extracted channels
output_dir = "./2025-05-09/extracted_channels/"
os.makedirs(output_dir, exist_ok=True)  # Create the directory if it doesn't exist

# Path to the multi-channel WAV file
for file in audio_files:
    file_path = os.path.join(DIR, file)

    angle_name = file.split('.')[0]
    if do_print:
        print(f"Processing file: {angle_name}")

    # Read the multi-channel WAV file
    audio_data, sample_rate = sf.read(DIR + file)

    # Check the shape of the audio data
    if do_print:
        print(f"Audio data shape: {audio_data.shape}")  # (samples, channels)

    # Extract individual channels
    num_channels = audio_data.shape[1]  # Number of channels
    channels = [audio_data[:, i] for i in range(num_channels)]

    # Save each channel as a separate WAV file
    for i, channel_data in enumerate(channels):
        output_file = os.path.join(output_dir + "channel_separation", angle_name+f"_{i + 1}.wav")  # Path to the output file
        sf.write(output_file, channel_data, sample_rate)
        if do_print:
            print(f"Saved channel {i + 1} to {output_file}")

2. Sorting of the angular extracted files into the separate channels

# List all extracted channel files separated by channel number
from natsort import natsorted

# Directory containing the extracted channels
extracted_channels_dir = "./2025-05-09/extracted_channels/channel_separation/"

# List all extracted channel files
channel_files = os.listdir(extracted_channels_dir)

# Filter out directories and non-audio files, keep only valid files
channel_files = [f for f in channel_files if os.path.isfile(os.path.join(extracted_channels_dir, f)) and f.endswith('.wav')]

# Sort the files naturally by the last part of their names (e.g., channel number)
sorted_channel_files = natsorted(channel_files, key=lambda x: int(x.split('_')[-1].split('.')[0]))

# Group files by the last part of their name (channel number)
grouped_files = {}

for file in sorted_channel_files:
    # Extract the channel number from the file name (e.g., "350_1.wav" -> "1")
    channel_number = int(file.split('_')[-1].split('.')[0])

    # Group files by channel number
    if channel_number not in grouped_files:
        grouped_files[channel_number] = []
    grouped_files[channel_number].append(file)

for i in range(len(grouped_files)):
    grouped_files[i+1].sort()

# Print grouped files
if do_print:
    for channel_number, files in grouped_files.items():
        print(f"Channel {channel_number}:")
        for f in files:
            print(f"  {f}")

3. Matched filtering:

• Every file is matched with the output
• First sweep is extracted and displayed
• First sweep is saved to a new directory to be analyzed
• avarage RMS of all the sweep in every recording is saved

import scipy
import scipy.stats

gain = 20  # dB

db_to_linear = lambda X: 10**(X/20)
dB = lambda X: 20 * np.log10(abs(np.array(X).flatten()))

extracted_channels_dir = "./2025-05-09/extracted_channels/channel_separation/"
saving_dir = "./2025-05-09/extracted_channels/"

# Define the matched filter function
def matched_filter(recording, chirp_template):
    filtered_output = np.roll(signal.correlate(recording, chirp_template, 'same', method='direct'), -len(chirp_template)//2)
    filtered_output *= signal.windows.tukey(filtered_output.size, 0.1)
    filtered_envelope = np.abs(signal.hilbert(filtered_output))
    return filtered_envelope

# Detect peaks in the matched filter output
def detect_peaks(filtered_output):
    peaks, properties = signal.find_peaks(filtered_output, prominence=0.1, distance=0.2 * sample_rate)
    return peaks

# Dictionary to store RMS values for all files
rms_values_dict = {}

channel_number = 0
for i in range(len(grouped_files)):
    channel_number += 1
    files = grouped_files[i+1]
    if do_print:
        print(f"\nProcessing Channel {channel_number}:\n")
    
        # Create a new figure for each channel
        fig, ax = plt.subplots(figsize=(15, 5))
        ax.set_title(f"Channel {channel_number}", fontsize=20)
        ax.set_xlabel("Seconds")
        ax.set_ylabel("Amplitude")
        ax.grid(True)

    for file in files[0:len(files)]:
        file_path = os.path.join(extracted_channels_dir, file)
        recording, sample_rate = sf.read(file_path)

        angle_name = int(file.split('_')[0])

        if 100 <= angle_name <= 260:  
            recording *= db_to_linear(-gain)
        rms_values = []

        # Apply matched filtering for each chirp template
        for idx_chirp, chirp_template in enumerate(chirp):
            filtered_output = matched_filter(recording, chirp_template)

            # Detect peaks
            peaks = detect_peaks(filtered_output)
            if do_print:
                print(f"Peaks detected in {file}: {len(peaks)} with {idx_chirp+1} chirp template")

                ax.plot(np.linspace(0, len(filtered_output) / sample_rate, len(filtered_output)), filtered_output, label=f"Filtered Output - {file}")
                #ax.plot(np.linspace(0, len(filtered_envelope) / sample_rate, len(filtered_envelope)), filtered_envelope, label=f"Filtered env - {file}", color='orange')
                ax.plot(np.linspace(0, len(recording) / sample_rate, len(recording)), recording, label=f"rec - {file}")
                #ax.legend(loc='upper right')
                #ax.set_xlim(5.15, 5.18)
                ax.set_xlim(2, 7)

                # Plot peaks on the filtered output
                ax.plot(peaks / sample_rate, filtered_output[peaks], "x", label="Detected Peaks")
                #ax.legend(loc='upper right')
            
            if len(peaks) > 0:
                # Extract the first sweep
                sweep_start = peaks[idx_chirp]
                sweep_end = sweep_start + len(chirp_template)
                sweep = recording[sweep_start:sweep_end]

                # Calculate RMS value of the first sweep
                rms = np.sqrt(np.mean(sweep**2))


                if do_print:
                    print(f"RMS value of sweep at peak {peaks[idx_chirp]/sample_rate} sec in {file}: {rms:.5f}")
                # Process each channel
                DIR_sweep = saving_dir + f"sweep_{idx_chirp+1}/"  # Directory to save the first sweeps
                os.makedirs(DIR_sweep, exist_ok=True)  # Create the directory if it doesn't exist

                sf.write(DIR_sweep + file, sweep, int(sample_rate))
                
                # Store RMS value in the dictionary
            else:
                if do_print:
                    print(f"No sweeps detected in {file} - Channel {channel_number}")
            # # Calculate RMS values for all detected peaks
                # if do_print:
                #     print(f"Average RMS value of all sweeps in {file}: {average_rms:.5f}")
            
            rms_values.append(rms)
            if do_print:
                print(f"rms_values in {file}: {rms_values}")
            

        rms_values_dict[file] = np.mean(rms_values)

# Print the dictionary of RMS values
if do_print:
    print("\nRMS Values for All Files:")
    for file, rms_value in rms_values_dict.items():
        print(f"{file}: {rms_value:.5f}")



plt.show(block = False)

4. Filtered sweeps in the files can be visulised

# Plot all angles
angles = [file.split('_')[0] for file in files]  # Extract angle names from filenames

sweep_to_plot = 5 # Sweep number to plot
DIR_sweep = saving_dir + f"sweep_{sweep_to_plot}/"  # Directory to save the first sweeps

channel_number = 0
for i in range(len(grouped_files)):
    channel_number += 1
    files = grouped_files[i+1]
    # Plot all angles, skipping '360'
    fig1, axs = plt.subplots(9, 4, figsize=(15, 30), sharey=True)
    angles = [file.split('_')[0] for file in files]  # Extract angle names from filenames

    idx_to_plot = 0
    for idx, file in enumerate(files):
        if angles[idx] == '360':
            continue  # Skip the 360 angle

        file_path = os.path.join(DIR_sweep, file)
        audio, sample_rate = sf.read(file_path)

        rms = np.sqrt(np.mean(audio**2))
        rms_db = 20 * np.log10(rms)
        
        row = idx_to_plot // 4
        col = idx_to_plot % 4
        
        ax = axs[row, col]
        ax.plot(np.linspace(0, len(audio) / sample_rate, len(audio)), audio)
        ax.set_title(f"Angle: {angles[idx]} degrees\nfile: {file}")  # Use extracted angle name with units
        ax.set_xlabel("Time (s)")
        ax.set_ylabel(f"Amplitude")
        ax.grid(True)
        ax.legend([f'RMS: {rms:.5f}\nRMS: {rms_db:.5f} dB'], loc='upper left')

        idx_to_plot += 1
    # Add suptitle after plotting all subplots
    plt.suptitle(f"Channel {channel_number}: Sweep number {sweep_to_plot}; Directory: {DIR_sweep}", fontsize=20, y=1)
    plt.tight_layout()  # Adjust layout to make room for suptitle
    plt.show(block = False)

    ax.legend()

plt.show(block = False)

# Print the dictionary of RMS values
print("\nRMS Values for All Files:")
for file, rms_value in rms_values_dict.items():
    print(f"{file}: {rms_value:.5f}")

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

RMS Values for All Files:
000_1.wav: 0.16406
010_1.wav: 0.17057
020_1.wav: 0.18191
030_1.wav: 0.19529
040_1.wav: 0.17347
050_1.wav: 0.15782
060_1.wav: 0.13673
070_1.wav: 0.11877
080_1.wav: 0.10848
090_1.wav: 0.09653
100_1.wav: 0.01894
110_1.wav: 0.01753
120_1.wav: 0.01429
130_1.wav: 0.01405
140_1.wav: 0.01331
150_1.wav: 0.01257
160_1.wav: 0.01323
170_1.wav: 0.01142
180_1.wav: 0.00992
190_1.wav: 0.01139
200_1.wav: 0.01545
210_1.wav: 0.01743
220_1.wav: 0.01819
230_1.wav: 0.01828
240_1.wav: 0.01828
250_1.wav: 0.01965
260_1.wav: 0.02159
270_1.wav: 0.10389
280_1.wav: 0.11129
290_1.wav: 0.11921
300_1.wav: 0.12378
310_1.wav: 0.12887
320_1.wav: 0.13464
330_1.wav: 0.14021
340_1.wav: 0.14579
350_1.wav: 0.15736
360_1.wav: 0.16406
000_2.wav: 0.20325
010_2.wav: 0.21812
020_2.wav: 0.21334
030_2.wav: 0.18878
040_2.wav: 0.16972
050_2.wav: 0.15631
060_2.wav: 0.13469
070_2.wav: 0.11697
080_2.wav: 0.10072
090_2.wav: 0.09700
100_2.wav: 0.01981
110_2.wav: 0.01952
120_2.wav: 0.01704
130_2.wav: 0.01365
140_2.wav: 0.01250
150_2.wav: 0.01236
160_2.wav: 0.01221
170_2.wav: 0.01253
180_2.wav: 0.01119
190_2.wav: 0.01047
200_2.wav: 0.01047
210_2.wav: 0.01461
220_2.wav: 0.01861
230_2.wav: 0.01956
240_2.wav: 0.01933
250_2.wav: 0.01885
260_2.wav: 0.01989
270_2.wav: 0.10592
280_2.wav: 0.11867
290_2.wav: 0.13328
300_2.wav: 0.14285
310_2.wav: 0.15088
320_2.wav: 0.15521
330_2.wav: 0.16166
340_2.wav: 0.17123
350_2.wav: 0.19221
360_2.wav: 0.20325
000_3.wav: 0.19741
010_3.wav: 0.18972
020_3.wav: 0.18283
030_3.wav: 0.17164
040_3.wav: 0.16237
050_3.wav: 0.15420
060_3.wav: 0.13550
070_3.wav: 0.11564
080_3.wav: 0.09741
090_3.wav: 0.08967
100_3.wav: 0.01847
110_3.wav: 0.01892
120_3.wav: 0.01827
130_3.wav: 0.01633
140_3.wav: 0.01243
150_3.wav: 0.00996
160_3.wav: 0.01046
170_3.wav: 0.01074
180_3.wav: 0.01141
190_3.wav: 0.01058
200_3.wav: 0.01033
210_3.wav: 0.00975
220_3.wav: 0.01110
230_3.wav: 0.01532
240_3.wav: 0.01716
250_3.wav: 0.01762
260_3.wav: 0.01735
270_3.wav: 0.09088
280_3.wav: 0.10333
290_3.wav: 0.11953
300_3.wav: 0.13304
310_3.wav: 0.14351
320_3.wav: 0.15778
330_3.wav: 0.17248
340_3.wav: 0.18480
350_3.wav: 0.19983
360_3.wav: 0.19741
000_4.wav: 0.18865
010_4.wav: 0.17766
020_4.wav: 0.17070
030_4.wav: 0.17119
040_4.wav: 0.17034
050_4.wav: 0.16293
060_4.wav: 0.15045
070_4.wav: 0.13402
080_4.wav: 0.11502
090_4.wav: 0.10315
100_4.wav: 0.01898
110_4.wav: 0.01831
120_4.wav: 0.01815
130_4.wav: 0.01796
140_4.wav: 0.01735
150_4.wav: 0.01517
160_4.wav: 0.01122
170_4.wav: 0.01074
180_4.wav: 0.01097
190_4.wav: 0.01126
200_4.wav: 0.01135
210_4.wav: 0.01118
220_4.wav: 0.01167
230_4.wav: 0.01205
240_4.wav: 0.01463
250_4.wav: 0.01670
260_4.wav: 0.01747
270_4.wav: 0.09472
280_4.wav: 0.10363
290_4.wav: 0.11743
300_4.wav: 0.13053
310_4.wav: 0.14038
320_4.wav: 0.15881
330_4.wav: 0.18626
340_4.wav: 0.19816
350_4.wav: 0.19454
360_4.wav: 0.18865
000_5.wav: 0.16082
010_5.wav: 0.15824
020_5.wav: 0.15822
030_5.wav: 0.15582
040_5.wav: 0.15112
050_5.wav: 0.14293
060_5.wav: 0.13698
070_5.wav: 0.12922
080_5.wav: 0.11938
090_5.wav: 0.11153
100_5.wav: 0.02093
110_5.wav: 0.01974
120_5.wav: 0.01872
130_5.wav: 0.01740
140_5.wav: 0.01672
150_5.wav: 0.01517
160_5.wav: 0.01476
170_5.wav: 0.01277
180_5.wav: 0.01090
190_5.wav: 0.01111
200_5.wav: 0.01230
210_5.wav: 0.01185
220_5.wav: 0.01165
230_5.wav: 0.01234
240_5.wav: 0.01353
250_5.wav: 0.01525
260_5.wav: 0.01683
270_5.wav: 0.10021
280_5.wav: 0.11238
290_5.wav: 0.12467
300_5.wav: 0.13638
310_5.wav: 0.14822
320_5.wav: 0.17008
330_5.wav: 0.17245
340_5.wav: 0.16342
350_5.wav: 0.17119
360_5.wav: 0.16082

5. RMS values of the overall recording for each microphone of the array are calculated and displayed

# RMS values of the overall recording for each channel and each angle

num_channels = len(grouped_files)
fig_polar, axs_polar = plt.subplots(num_channels, 1, figsize=( 5,25), subplot_kw={'projection': 'polar'})
fig_polar.suptitle("Overall RMS Values", fontsize=16)

fig_linear, ax_linear = plt.subplots(figsize=(10, 10))
fig_linear.suptitle("Overall RMS Values for All Channels", fontsize=16)

for i in range(num_channels):
    channel_number = i + 1
    files = grouped_files[channel_number]
    
    rms_values = []
    rms_values_norm_db = []
    angles = []
    
    for file in files:

        rms = rms_values_dict[file]
        
        rms_values.append(rms)

        rms_values_norm = rms_values / rms_values[0]
        rms_values_norm_db = 20 * np.log10(rms_values_norm)

        angle_name = file.split('_')[0]
        angles.append(int(angle_name))
    
    # Convert angles to radians
    angles_rad = np.radians(angles)
    
    # Plot RMS values in polar plot
    ax_polar = axs_polar[i] if num_channels > 1 else axs_polar
    ax_polar.plot(angles_rad, rms_values_norm_db, linestyle='-')
    ax_polar.set_title(f"Mic {channel_number}")
    ax_polar.set_theta_zero_location("N")  # Set 0 degrees to North
    ax_polar.set_theta_direction(-1)  # Set clockwise direction
    ax_polar.set_xticks(np.linspace(0, 2 * np.pi, 18, endpoint=False))  # Set angle ticks
    ax_polar.set_xlabel("Angle (degrees)") 
    ax_polar.tick_params(axis='y', labelsize=8)
    ax_polar.set_ylabel("RMS Value dB", position=(0, 1), ha='left')
    ax_polar.set_rlabel_position(0)
    ax_polar.set_yticks(np.arange(-30, 3, 2))  # Set y-ticks

    # Plot RMS values in linear plot
    ax_linear.plot(angles, rms_values_norm_db, marker='.', linestyle='-', label=f"Mic {channel_number}")
    ax_linear.set_xlabel("Angle (degrees)")
    ax_linear.set_xticks(np.linspace(0, 380, 19, endpoint=False))  # Set angle ticks
    ax_linear.set_ylabel("RMS Value dB")
    ax_linear.set_yticks(np.arange(-30, 5, 1))  # Set y-ticks
    ax_linear.legend()
    ax_linear.grid(True)

fig_polar.tight_layout()
plt.show(block = False)

6. Frequency analysis of every channel up to 40 KHz

import soundfile as sf
from scipy import fft

sweep_to_plot = 1  # Microphone number to plot
DIR = f"./2025-05-09/extracted_channels/sweep_{sweep_to_plot}/"  # Directory containing the audio files
# Central frequencies of the bands
central_freq = np.array([4e3, 6e3, 8e3, 10e3, 12e3, 14e3, 16e3, 18e3, 20e3, 22e3])
BW = 1e3  # Bandwidth of the bands
linestyles = ["-", "--", "-.", ":", "-", "--"]  # Line styles for the plot
# Group central frequencies into 3 sets of 6 bands each
num_bands_per_plot = 6
central_freq_sets = [central_freq[i * num_bands_per_plot:(i + 1) * num_bands_per_plot] for i in range(3)]
# Number of microphones
num_mics = num_channels
# Plot for each microphone
for mic in range(1, num_mics + 1):
    files = grouped_files[mic]
    angles = [int(file.split('_')[0]) for file in files]  # Extract angles from filenames
    # Create a figure with 3 polar subplots
    fig, axes = plt.subplots(1, 2, subplot_kw={"projection": "polar"}, figsize=(10, 6))
    plt.suptitle(f"Polar Frequency Response, sweep {sweep_to_plot}- Microphone {mic}", fontsize=20)
    for ax_idx, ax in enumerate(axes):
        ii = 0
        for fc in central_freq_sets[ax_idx]:
        
            audio_patt = []
            for file in files:
                file_path = os.path.join(DIR, file)
                audio, fs = sf.read(file_path)
                # Compute FFT
                audio_freq = fft.fft(audio, n=2048)
                audio_freq = audio_freq[:1024]
                freqs = fft.fftfreq(2048, 1 / fs)[:1024]
                # Compute mean radiance in the band
                band_mean = np.mean(np.abs(audio_freq[(freqs > fc - BW) & (freqs < fc + BW)]))
                audio_patt.append(band_mean)
            # Normalize and plot
            audio_patt_norm = audio_patt / audio_patt[0] # Normalize the radiance
            audio_patt_norm_dB = 20 * np.log10(audio_patt_norm) # Convert the radiance to dB
        
            if fc >= 10e3:
                label = f"{fc / 1e3:.0f} kHz"
            else:
                label = f"{fc / 1e3:.0f} kHz"
        
            ax.plot(np.deg2rad(angles), audio_patt_norm_dB, label=label, linestyle=linestyles[ii])
            ii +=1
        # Configure polar plot
        ax.legend(loc="upper right")
        ax.set_theta_offset(np.pi / 2)
        ax.set_theta_zero_location("N")  # Set 0 degrees to North
        ax.set_theta_direction(-1)  # Set clockwise direction
        ax.set_xticks(np.linspace(0, 2 * np.pi, 18, endpoint=False))  # Set angle ticks
        ax.set_yticks(np.linspace(-40, 0, 7))
        ax.set_xlabel("Angle (degrees)")
        ax.set_ylabel("RMS Value dB", position=(0, 1), ha='left')
        ax.set_rlabel_position(0)
    plt.tight_layout()
    plt.show()

7. Noise from every recorded file is trimmed and saved

channel_number = 0
for i in range(len(grouped_files)):
    channel_number += 1
    files = grouped_files[i+1]
    if do_print:
        print(f"\nProcessing Channel {channel_number}:\n")
    
    for file in files[0:len(files)]:
        file_path = os.path.join(extracted_channels_dir, file)
        recording, sample_rate = sf.read(file_path)

        if 90 <= angle_name <= 250:  
            recording *= db_to_linear(-gain)

        angle_name = int(file.split('_')[0])

        start = int(sample_rate * 0.5)  # 100ms
        end = start + int(sample_rate * 0.005)  # 100ms
        noise = recording[start:end]

        DIR_noise = saving_dir + "ground_noises/"  # Directory to save the first sweeps
        os.makedirs(DIR_noise, exist_ok=True)  # Create the directory if it doesn't exist

        sf.write(DIR_noise + file, noise, int(sample_rate))

8. Trimmed noises can be visualized in along 360 degrees for every channel

# Plot all angles, skipping '360'
angles = [file.split('_')[0] for file in files]  # Extract angle names from filenames

sweep_to_plot = 5

channel_number = 0
for i in range(len(grouped_files)):
    channel_number += 1
    files = grouped_files[i+1]
    # Plot all angles, skipping '360'
    fig1, axs = plt.subplots(9, 4, figsize=(15, 30), sharey=False)
    angles = [file.split('_')[0] for file in files]  # Extract angle names from filenames

    idx_to_plot = 0
    for idx, file in enumerate(files):
        if angles[idx] == '360':
            continue  # Skip the 360 angle

        file_path = os.path.join(DIR_noise, file)
        noise, sample_rate = sf.read(file_path)

        rms = np.sqrt(np.mean(noise**2))
        rms_db = 20 * np.log10(rms)
        
        row = idx_to_plot // 4
        col = idx_to_plot % 4
        
        ax = axs[row, col]
        ax.plot(np.linspace(0, len(noise) / sample_rate, len(noise)), noise)
        ax.set_title(f"Angle: {angles[idx]} degrees\nfile: {file}")  # Use extracted angle name with units
        ax.set_xlabel("Time (s)")
        ax.set_ylabel(f"Amplitude")
        ax.grid(True)
        ax.legend([f'RMS: {rms:.7f}\nRMS: {rms_db:.5f} dB'], loc='upper left')

        idx_to_plot += 1
    # Add suptitle after plotting all subplots
    plt.suptitle(f"Channel {channel_number}; Directory: {DIR_noise}", fontsize=20, y=1)
    plt.tight_layout()  # Adjust layout to make room for suptitle
    plt.show(block = False)

    ax.legend()

plt.show(block = False)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

9. Signal to Noise Ratio is calculated and displayed for all the files, considering the +20 dB appllied from 100 to 260 degrees

#%% overall SNR computation

def calculate_snr(audio, noise):

    # Compute signal power (RMS value squared)
    P_signal = np.mean(audio**2)

    # If noise segment is provided, extract noise power
    P_noise = np.mean(noise**2)

    # Compute SNR in dB
    SNR = 10 * np.log10(P_signal / P_noise)

    # Noise floor in dBFS (full scale)
    noise_floor = 10 * np.log10(P_noise)

    return SNR, noise_floor

sweep_to_plot = 5  # Sweep used for SNR calculation

DIR = f"./2025-05-09/extracted_channels/sweep_{sweep_to_plot}/"  # Directory containing the audio files
fig_linear, ax_linear = plt.subplots(figsize=(10, 5))
fig_linear.suptitle("Overall RMS Values for All Channels", fontsize=16)

SNR_values_dict = {}
num_mics = num_channels
# Plot for each microphone
for mic in range(1, num_mics + 1):
    files = grouped_files[mic]
    angles = [int(file.split('_')[0]) for file in files]  # Extract angles from filenames
    # Create a figure with 3 polar subplots

    plt.suptitle(f"SNR using sweep {sweep_to_plot}", fontsize=20)

    for file in files:
        file_path = os.path.join(DIR, file)
        audio, fs = sf.read(file_path)

        noise_file_path = os.path.join(DIR_noise, file)
        noise, fs = sf.read(noise_file_path)

        # Calculate SNR
        SNR, noise_floor = calculate_snr(audio, noise)
        if do_print:
            print(f"SNR for {file}: {SNR:.2f} dB, Noise Floor: {noise_floor:.2f} dBFS")
        SNR_values_dict[file] = SNR
        # Normalize and plot
        # audio_patt_norm = SNR / SNR[0] # Normalize the radiance
        SNR_dB = 20 * np.log10(SNR) # Convert the radiance to dB


    # Extract SNR values for the current mic in the order of angles/files
    snr_values = [SNR_values_dict[file] for file in files]

    # Plot SNR values in linear plot
    ax_linear.plot(angles, snr_values, marker='.', linestyle='-', label=f"Mic {mic}")
    ax_linear.set_xlabel("Angle (degrees)")
    ax_linear.set_xticks(np.linspace(0, 380, 19, endpoint=False))  # Set angle ticks
    ax_linear.set_ylabel("SNR (dB)")
    ax_linear.set_yticks(np.arange(0, 16, 2))  # Adjust y-ticks for SNR range
    ax_linear.legend()
    ax_linear.grid(True)

10. Sensitivity

import numpy as np 
import soundfile as sf
import matplotlib.pyplot as plt 
from utilities import *
import scipy.signal as sig

# Load the substitution-calibration audio (calibration and target mic)
fs_gras = sf.info('./2025-05-09/ref_tone_gras2025-05-09_11-59-29.wav').samplerate
gras_pbk_audio, fs_gras = sf.read('./2025-05-09/ref_tone_gras2025-05-09_11-59-29.wav',
                         start=int(fs_gras*3.530), stop=int(fs_gras*3.533))

fs_SPH0645 = sf.info('./2025-05-09/extracted_channels/channel_separation/000_3.wav').samplerate
SPH0645_pbk_audio, fs_SPH0645 = sf.read('./2025-05-09/extracted_channels/channel_separation/000_3.wav',
                               start=int(fs_SPH0645*4.334),  stop=int(fs_SPH0645*4.337))
# resampling 
gras_pbk_audio = sig.resample(gras_pbk_audio, int(fs_SPH0645*0.003))


plt.figure(figsize=(10,5))
plt.plot(np.linspace(0,len(gras_pbk_audio)/fs_SPH0645, len(gras_pbk_audio)) ,gras_pbk_audio)
plt.title('GRAS playback signal')
plt.figure(figsize=(10,5))
plt.plot(np.linspace(0,len(SPH0645_pbk_audio)/fs_SPH0645, len(SPH0645_pbk_audio)), SPH0645_pbk_audio)
plt.title('SPH0645 playback signal')

# Load the 1 Pa reference tone 
gras_1Pa_tone, fs_gras = sf.read('./2025-05-09/gras_1Pa_ch9_30dB_chA_20dB.wav', start=int(fs_gras*0.5),
                        stop=int(fs_gras*1.5))


gras_1Pa_tone = sig.resample(gras_1Pa_tone, fs_SPH0645*1)

plt.figure(figsize=(15,5))
plt.plot(np.linspace(0,len(gras_1Pa_tone)/fs_SPH0645, len(gras_1Pa_tone)), gras_1Pa_tone)
plt.title('SPH0645 playback signal')

Text(0.5, 1.0, 'SPH0645 playback signal')

# Calibration mic: Calculate the rms_Pascal of the 1 Pa calibration tone
rms_1Pa_tone = rms(gras_1Pa_tone)
print(f'The calibration mic has a sensitivity of {np.round(rms_1Pa_tone,3)}rms/Pa. RMS relevant only for this ADC!')

# Now measure mic RMS over all frequency bands
gras_centrefreqs, gras_freqrms = calc_native_freqwise_rms(gras_pbk_audio, fs_SPH0645)

# Convert from RMS to Pascals (rms equivalent) since we know the GRAS sensitivity
gras_freqParms = gras_freqrms/rms_1Pa_tone # now the levels of each freq band in Pa_rms

plt.figure()
a0 = plt.subplot(211)
plt.plot(gras_centrefreqs, gras_freqParms)
plt.ylabel('Pressure_rmseqv., Pa', fontsize=12)
plt.title('GRAS mic recording of playback')
plt.subplot(212, sharex=a0)
plt.plot(gras_centrefreqs, pascal_to_dbspl(gras_freqParms))
plt.xlabel('Frequencies, Hz', fontsize=12);
plt.ylabel('Sound pressure level,\n dBrms SPL re 20$\mu$Pa', fontsize=12)

# Target microphone. Here we'll cover the case where we only get an RMS/Pa
# sensitivity. The other option is to calculate a mV/Pa sensitivity - which allows
# you to use the mic aross different ADCs - but also needs more info on the ADC specs

SPH0645_centrefreqs, SPH0645_freqrms = calc_native_freqwise_rms(SPH0645_pbk_audio, fs_SPH0645)
plt.figure()

a0 = plt.subplot(211)
plt.plot(SPH0645_centrefreqs, SPH0645_freqrms)
plt.ylabel('a.u. rmseqv.', fontsize=12)
plt.title('SPH0645 mic recording of playback')
plt.subplot(212, sharex=a0)
plt.plot(SPH0645_centrefreqs, dB(SPH0645_freqrms))
plt.xlabel('Frequencies, Hz', fontsize=12);
plt.ylabel('dBrms a.u.', fontsize=12)


# Now let's calculate the RMS/Pa sensitivity using the knowledge from the 
# calibration mic
SPH0645_sensitivity = np.array(SPH0645_freqrms)/np.array(gras_freqParms)

plt.figure()
a0 = plt.subplot(211)
plt.plot(SPH0645_centrefreqs, SPH0645_sensitivity)
plt.ylabel('a.u. RMS/Pa', fontsize=12)
plt.title('Target mic sensitivity')
plt.subplot(212, sharex=a0)
plt.plot(SPH0645_centrefreqs, dB(SPH0645_sensitivity))
plt.xlabel('Frequencies, Hz', fontsize=12);
plt.ylabel('dB a.u. rms/Pa', fontsize=12)
# plt.ylim(-60,0)

The calibration mic has a sensitivity of 0.043rms/Pa. RMS relevant only for this ADC!

Text(0, 0.5, 'dB a.u. rms/Pa')

# We now have the target mic sensitivity - how do we use it to calculate the
# actual dB SPL? 

# Here we load a separate 'recorded sound' - a 'validation' audio clip let's call it 

recorded_sound, fs_SPH0645 = sf.read('./2025-05-09/extracted_channels/channel_separation/000_3.wav',
                               start=int(fs_SPH0645*4.537),  stop=int(fs_SPH0645*4.541))

# Also load the 'validation' calibration mic recording of the same sound
gras_rec, fs_gras = sf.read('./2025-05-09/ref_tone_gras2025-05-09_11-59-29.wav',
                         start=int(fs_gras*3.733), stop=int(fs_gras*3.737))
gras_rec = sig.resample(gras_rec, int(fs_SPH0645*0.004))

plt.figure(figsize=(10,5))
plt.plot(np.linspace(0,len(gras_rec)/fs_SPH0645, len(gras_rec)) ,gras_rec)
plt.title('GRAS playback signal')
plt.figure(figsize=(10,5))
plt.plot(np.linspace(0,len(recorded_sound)/fs_SPH0645, len(recorded_sound)), recorded_sound)
plt.title('SPH0645 playback signal')

Text(0.5, 1.0, 'SPH0645 playback signal')

# And finally let's check that the Sennheiser calibration makes sense
# using a sound that we didn't use to calculate the sensitivity
# If the length of the recorded target mic audio here is not the same as the calibration audio. 
#  then you'll need to interpolate the microphone sensitivity using interpolate_freq_response in the
# utilities.py module
recsound_centrefreqs, freqwise_rms = calc_native_freqwise_rms(recorded_sound, fs_SPH0645)
interp_sensitivity = interpolate_freq_response([SPH0645_centrefreqs, SPH0645_sensitivity],
                          recsound_centrefreqs)
freqwise_Parms = freqwise_rms/interp_sensitivity # go from rms to Pa(rmseq.)
freqwiese_dbspl = pascal_to_dbspl(freqwise_Parms)


gras_centrefreqs, gras_freqrms = calc_native_freqwise_rms(gras_rec, fs_SPH0645)
gras_Pa = gras_freqrms/rms_1Pa_tone
gras_dbspl = pascal_to_dbspl(gras_Pa)

plt.figure()
plt.plot(gras_centrefreqs,gras_dbspl, label='gras')
plt.plot(recsound_centrefreqs,freqwiese_dbspl, label='SPH0645')
plt.ylabel('dBrms SPL, re 20$\mu$Pa', fontsize=12)
plt.xlabel('Frequency, Hz', fontsize=12)
plt.legend()
plt.show()

# Now we know the sensitivity of the target mic - let's finally calculate
# the dB SPL of the recorded sound!
# We rely on combining the Pa rms of all relevant frequencies 
# e.g. see https://electronics.stackexchange.com/questions/642109/summing-rms-values-proof

frequency_band = [0.2e3, 15e3] # min, max frequency to do the compensation Hz

# Choose to calculate the dB SPL only for the frequency range of interest.
# Target mic
tgtmic_relevant_freqs = np.logical_and(recsound_centrefreqs>=frequency_band[0],
                                recsound_centrefreqs<=frequency_band[1])
total_rms_freqwise_Parms = np.sqrt(np.sum(freqwise_Parms[tgtmic_relevant_freqs]**2))

# Ground truth GRAS mic audio of the same sound. Here we use only the relevant 
# frequency band of the recorded sweep
gras_relevant_freqs = np.logical_and(gras_centrefreqs>=frequency_band[0],
                                gras_centrefreqs<=frequency_band[1])

#%% This is one way to do it - use the RAW audio, and calculate its rms
# since we the overall flat sensitivity of the GRAS
gras_overallaudio_Parms = rms(gras_rec)/rms_1Pa_tone
# This is the other way to do it by combining RMS values - like we did for the Sennheiser
gras_totalrms_Parms = np.sqrt(np.sum(gras_Pa[gras_relevant_freqs]**2))

print(f'GRAS dBrms SPL measures:{pascal_to_dbspl(gras_overallaudio_Parms)}, {pascal_to_dbspl(gras_totalrms_Parms)}')
print(f'SPH0645 dBrms SPL measure: {pascal_to_dbspl(total_rms_freqwise_Parms)}')

GRAS dBrms SPL measures:[101.95870282], [101.6232441]
SPH0645 dBrms SPL measure: [101.62078077]

Final notes

point 4: the 5 sweeps ina single recording are individually extracted by correlating the correspondent template chirp. Each individual extracted chirp is saved and used later on for comparative analysis. Longer and shorter chirp do not change the overall results. RMS calculated values, however are still a mean of all the chirps. Sweep 5 is selected as example of the extracted signals from recordings

point 5: the overall RMS around the array on th robot hgas a drop in level at the back. Array was rotated around the central mic = 3, which is the only symmetrical in the rms resulting plot. However the asymmetry matches the fact that the other microphone were moving closer of further away from the sound source.

point 9: SNR has been kept bigger than 5 dB by uincreasing the level of the output from 100 to 260 degrees. The SNR however is not higher due to distortion, caused by clipping, seen from the mics if the output was further increased.

point 10: Mic sensitivity of SPH0645 MEMS from Adafriut is characterized. The experimental conditions match with the approximate 101 dBrms SPL obtained.